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A classical problem in elasticity theory involves an inhomogeneity embedded in a material of given 
stress and shear moduli. The inhomogeneity is a region of arbitrary shape whose stress and shear 
moduli differ from those of the surrounding medium. In this paper we present a new, semi-analytic 
method for finding the stress tensor for an infinite plate with such an inhomogeneity. The solution 
involves two conformal maps, one from the inside and the second from the outside of the unit circle 
to the inside, and respectively outside, of the inhomogeneity. The method provides a solution by 
matching the conformal maps on the boundary between the inhomogeneity and the surrounding 
material. This matching converges well only for relatively mild distortions of the unit circle due 
to reasons which will be discussed in the article. We provide a comparison of the present result to 
known previous results. 



I. INTRODUCTION 

Elasticity theory in homogeneous materials is a well 
developed subject. Much less is known about inhomoge- 
neous materials where the solution of the basic equations 
of elasticity becomes very involved. In this paper we fo- 
cus on a material which consists of one finite area (of 
arbitrary shape) in which a material with given elastic 
properties is embedded in an infinite sheet of material 
of different elastic properties. This situation is known 
as "elastic inhomogeneity" and it appears in a variety of 
solid mechanical contexts. Previous studies have concen- 
trated on solving this problem for the relatively symmet- 
ric case of an ellipse jl| where it was solved analytically. 
In other cases the problem was solved for small pertur- 
bations to the circle @. Mathematically the problem is 
set as follows, and see Fig. [TJ A patch of material of 
type 1 occupies an area fl and is delineated by a sharp 
boundary which will be denoted dfl. The rest of the infi- 
nite plane is made of material of type 2. The material is 
subjected to forces at infinity (and see below the precise 
boundary conditions), and is therefore deformed. Before 
the deformation each point of the material is assigned a 
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FIG. 1: The region Q. 



point r in the two-dimensional plane. The forces at in- 
finity result in a displacement of the material points to a 
new equilibrium position r' . The displacement field u(r) 
is defined as 

u(r) = r' — r . (1) 

The strain field is defined accordingly as 

eij ee [diUj + djUi]/2 . (2) 

In the context of linear elasticity in isotropic materials 
one then introduces the stress field according to Hook's 
law 

°y = 2/1^6^ + XiS^ekk (3) 

Where Ai = MiAoirr — !)■ Mi an d v% take on different 
values fii, v\ in Q and /i2, V2 in the rest of the material. 
In equilibrium the stress tensor should be divergencclcss 
= at each point in the sheet. By defining the stress 
(or Airy) potential U: 

d 2 U _ d 2 U _ d 2 U 

Uxx ~ ~dt ' ,(Txv ~~dxdy~ ' 1<7yy ~ fe 2 " ' ( ) 

the former equation for the stress tensor becomes a par- 
tial differential equation for the stress potential: 

V 2 V 2 U{x,y) = . (5) 

This equation, which is known as the bi-Laplace or the 
bi-harmonic equation is conveniently solved as a non- 
analytic combination of analytic funcitons. To this aim 
we introduce the complex notation z ee x + iy, and note 
the general solutions of Eq. ^ in the form 

U(x,y)=M[z<p(z)+ri(z)], (6) 

where z ee x — iy and <p(z) and rj(z) are any two holo- 
morphic functions. What remains to do in any particular 
problem is to find the unique analytic functions such that 
the stress tensor satisfies the boundary conditions. This 
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stress tensor is determined by the two holomorphic func- 
tions as: 



'yy 



(x,y) = K[2${z)+0'{z)+Tf'(z)\ 



(p'{z) + z(p"{ 
v xx {x,y) = ®.[2Cp'{z)-z(p"{z)- //'(:)] 
Oxy{x,y) = $S[z(p"(z) + r}"(z)]. 

We define for convenience: 



1p(z) = T)'(z) 



And then: 



(7) 



(8) 



on and Cp 2 and il>2 which are defined on C\f2. Never- 
theless we will demand continuity of the physical fields. 
In particular the normal force 



er • n = o~ x 



(11) 



and the displacement u(r) must be continues across the 
interface (by Newton's third law) in the absence of sur- 
face tension. Therefore the continuity conditions are: 



7 (D +io -(i) = CT (2) 

' xn 1 yn xn 



7 (2) 



(12) 



a y y(x,y) = $t[2tp'(z) + z(p"{z) + V>'(z)] 
a xx {x,y) = ?k[2(p , {z)-z(p"{z)-i>'{z)] 

v xy (x,y) = ^[z<p"(z) + i>'(z)]. 



(9) 



Note that the stress tensor is determined by derivatives 
of the holomorphic functions, and not by the functions 
themselves. This leaves us with some freedom, since the 
functions can be changed with the following gauge: 



& — > <Pi + iCiZ + a>i + if3, 

4>i — > 4>i + ii + iSi , i> 



(10) 



As we shall see below, not all these gauge freedoms are 
true freedoms once we introduce the boundary and conti- 
nuity conditions. Since the elastic properties are different 
inside and outside f2, the potential functions will be dif- 
ferent in the two regions: (p.\ and which are defined 



u^+iuV=uW+iuW . (13) 



y 



The continuity conditions for the stress, can be rewritten 
as: 



d_ fdlh , dUi 
ds \ dx dy 



£ [sr+i-sr =3i -ar (14) 



dU 2 



.dUo 



ds \ dx 



dy 



or, after integrating: 



dUi +i dUi _ d\h | .dU 2 | e 
dx dy dx dy 



(15) 



where C is a complex constant of integration. In terms of 
the holomorphic functions, the condition (fT2|) translates 
to H: 



(p {1) (z) + zif'W (z) + $W (z) = ^ (2) (z) + z(p'W (z) + V^ (2) (*) + C , 
and the condition fT3"l) becomes: 



(16) 



[ Kl ^ (z) - zcp'W (z) - i/>W (z)} [k 2 ^ (z) - ztp'V) (z) - ^(2) [z)\ 



Hi 



where Kj = (3 — + Vi)- 



In addition to these continuity conditions on d£l we need 
to specify boundary conditions at infinity. We choose 



& XX (oo) = , &yy{00) 



a xy (oo)=0. (18) 



/'2 



(17) 
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side and the point at infinity is outside fl. Since the 
stress functions are holomorphic in their domains of def- 
inition, we can expand them in the appropriate laurent 
series, which for the functions with superscript (1) is of 
the form 



II. SOLUTION IN TERMS OF CONFORMAL 
MAPS 

Solutions to the problem of finding the stress field out- 
side a given domain using conformal maps were described 
for example in [f|. Here we need to solve for the stress 
field both inside and outside the given domain. In the 
following we assume that the center of coordinates is in- 



^W(z) 



(19) 



i.e. we have no poles at the origin. For the outside do- 
main (functions with superscript (2)) the most general 
expansions in agreement with the boundary conditions 
(fT51) are of the form 
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and: 



,7,(2) 



~(2) , -(2) 



^)(z) = $W z + ^> +^(/ z + ^/ z * + ... (20) 

i.e we have a pole of order 1 at infinity. Accordingly, the 
leading terms of Eqs. (|20|) are determined by the bound- 
ary conditions. We now use one of the gauge freedoms 
to eliminate the imaginary part of (p^ and write: 



„7.( 2 ) 



(2) 



r,(2) . 



-(2) Coo r(2) CToo 



(21) 



The standard way to proceed [H would be to substi- 
tute the series expansions in the continuity conditions 
and find the linear equations that determine all the coeffi- 
cients by equating terms of the same order in z. However 
this cannot be done in general since the functions z n are 
not orthogonal on arbitrary contours dVL. To overcome 
this, one maps the regions f2 and C\0 into the interior 
and exterior of the unit circle, respectively. That is, we 
need two holomorphic, invertible (and thus conformal) 
functions, one is 



$(w) 



(22) 



which maps the exterior of the unit circle into C\f2, and 
the other is 



z = A(C) 



(23) 



which maps the unit disk into O. Since they are both 
invertible they have inverse functions which we denote 



and 



uj = $ _1 (z) 



(24) 



(25) 



Now we express the functions ip^ and ip^ in terms of 
uj and £ and then expand them on the boundary of the 
unit circle. This expansion will be a Fourier series where 
the powers of uj or £ satisfies the orthogonality relation: 



2tt 



(26) 



We have here used e 10 to represent either uj or £. The 
orthogonality allows us to equate the coefficients of the 
series term by term. We define: 



^{Z) = ^ (^- X {Z)) , ^(z)W = ^(2) ($-!(z)) 



(27) 



and: 



<pM{z) = <pV> (h-\z)) , 4> (1) {z) ee (A- 1 ^)) . 

We can expand tpi and ipi in terms of uj and £ on the unit 
circle. Since the original functions were holomorphic and 
meromorphic in the original domains, the functions: 

^ 2 \uj) = ^ 2) , V (2) M = ^ (2) ($M) . (29) 



^){Q = ^{K{Q), V (1) (C)^ (1) (A(C)) • (30) 

are holomorphic inside and outside the unit disc, respec- 
tively. Therefore we can expand in terms of uj and £: 



(31) 



^(2) ( w ) = ^C0 W + ^« + ^cj- 1 + ^cj- 2 + • • (32) 

We now assume that the map of the exterior domain, $, 
maps the point at infinity to infinity. That is, it will have 
a Laurent series on the form 



$(w) = Fiuj + F + F-iLU 1 + F-2U~ 



(33) 



From this we get the following relations (after substitut- 
ing and taking the limit uj — > oo): 



(2) _ p CToo ( 2 ) (Too 



(34) 



-.(») 



we can also use the last two freedoms to choose ip 

(2) (2) 

— Fo<p\ such that ip = 0. In the interior domain, the 
functions ip^ and tp^ also have five freedom. However, 
the requirement of continuity of the displacement field u 
across the boundary dQ removes three of these freedoms. 
This continuity was expressed by Eq. (|13[) . Applying 
the apparent gauge freedoms on the LHS of that equa- 
tion and then subtracting the resulting equation from Eq. 
(Ti"3"]) we find the three conditions 



Ci = , KlQll=7l, Kifli = -<5i 



(35) 



Using the remaining two freedoms, we can eliminate 
the constant term in the expansion of ip^ by setting 
■0Q 1 - 1 = 0. Note that this is possible only when we choose 
A(C) such that A(0) = 0. We may always define our 
mapping A such that this is satisfied. In terms of the 
conformal maps we transform the boundary conditions 
into 

(p<® (uj) + ^-^(p'Wjuj) + i/>W(u) , (36) 
$'(w) 



-1[ K1(/3 W(0- i^(i)( C )-v>W(C)] 

Mi A'(C) 

= - [K2<p {2) (w) - §^=V^~) - WH\ • (37) 
M2 $'(w) 
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III. METHOD OF SOLUTION 

At this point we need to substitute the expansions 
(|3T|l .([32 |) .([33 |l and an expansion similar to (|33| for A(£) 
into the equations ([55]) and (|3"Tj) and solve for the coef- 
ficients (p fi and ipk ^- To understand how to do this in 
principle we write the expanded equations (|3"Tj) and (|32|) 
in an abstract form 



E 



E 



(38) 



where pk are linear combinations of the coefficients ip 



(i) 



and tpn^ whereas q m are linear combinations of cpn J and 

("2") 

ipn . As this equation stands we cannot use the orthog- 
onality relation Eq. (|26[) . Therefore we expand moments 
of oj in terms of £ in the form 



(2) 



c(C) m = 



IS 



(39) 



We now insert this expression in Eq. (|38|) . 

oo oo oo 

E p^ k = E E ^ a »,-c n . (40) 



TTl= — oo n— — oo 



and equate the coefficients of same powers to achieve a 
set of linear equations for the coefficients of ip^ and 
The actual algebraic manipulations that are involved in 
reaching a finite set of linear equations are presented in 
the appendix. 

Needless to say, when we derive a finite set of equations 
we lose precision. To see this we note that to get the 
right number of equations for the number of unknowns 
(see Appendix) we need to truncate the summations on 
the LHS and the RHS of Eq. (@DJ at the same finite N, 
i.e. 



N 



N 



N 



E Pk ^ = E E ^m^mC" 



(41) 



-JV 



m=-N n=-N 



For a precisely circular inclusion this truncation intro- 
duces no loss of information. For this particular shape 
the expansion Eq. (|39[) has only one term with n = m, i.e. 
dn,m — 5 n .m- Obviously, when the inclusion shape devi- 
ates from the circle, the representation of ui m in Fourier 
space deviates from a delta function and it becomes more 
spread. An example of this phenomenon is presented in 
Figj2]for an inclusion in the form of an ellipse with aspect 
ratio of about 1.5. The upper panel shows the param- 
tcrization of the outer mapping as function of that of 
the inner, arg($ _1 (A(w))). In the lower panel we show 
the power spectrum |a„. m | 2 of the moments m — 1 and 
m = 25 of the function in the upper panel. If we trun- 
cate the expansion at the dashed line N = n = 25, we 
lose high frequency information for the higher moments. 
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FIG. 2: In the upper panel, we show the relation between the 
mapping of the inner and outer domains of the ellipse drawn 
in the inset. The ellipse is the same as the one used in the 
Figs. [5] and [7] Specifically, we have plotted the parametriza- 
tion of the outer mapping as function of that of the inner, 
arg(<E> (A(uj))). In the lower panel we show the power spec- 
trum |on, m | 2 of the moments m = 1 and m — 25 of the 
function in the upper panel. 



This loss of informaion will lead to stress field calcula- 
tions which are less accurate. 

To see the difficulty in a pictorial way we can consider 
the field lines of the conformal mappings for an inclu- 
sion that is elongated in shape, see for example Fig. [3] 
The external field lines concentrate at the convex parts of 
the inclusion, whereas the internal field lines concentrate 
on the concave parts. It becomes increasingly difficult 
to match field lines since they make a large discontin- 
ues jump when we go from the interior to the exterior 
domain. Similarly, for the ellipse in Fig. [21 when we in- 
crease the aspect ratio of the ellipse, the slope in the steep 
parts of arg($ _1 (A(o;))) become even larger, requiring 
higher order frequencies in our expansions. Eventually 
for large aspect ratios, our method will break down. 



FIG. 3: Field lines of the two conformal mappings to the 
interior and exterior domains, respectively. 



FIG. 5: a xx evaluated along the positive real axis (ellipse) 
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FIG. 4: 

IV. OBTAINING THE CONFORMAL MAPS 

In all the calculations we assumed that the conformal 
maps $(w) and A(£) are available. For arbitrary inclu- 
sion shapes this is far from obvious, and special methods 
are necessary to obtain these maps. An efficient method 
to obtain the conformal map from the exterior of the unit 
circle to the exterior of an arbitrary given shape had been 
discussed in great detail in [6| . In the present case we use 
a slightly different method namely the geodesic algorithm 
This method, like the former one, is based on the it- 
erations of a generic conformal map 7^ defined by a set 
of parameters £. We then construct the conformal map 
to an arbitrary shape by an appropriate choice of pa- 
rameters £. In the geodesic algorithm, we discretize the 
interface of the inclusion by a sequence of points {zk}^ =Q - 
The points appear sequentially in the positive direction 
of the interface. We now brieflysummarize how to con- 
struct the conformal map (see Q)- First we construct 
iteratively the inverse map that brings the interior of the 



inclusion to the upper half plane and the interface to the 
real axis. The conformal map to the shape then follows 
directly from the inverse. The construction is done in 
three steps. In the first step, we move one point to infin- 
ity and another to the center of coordinates, e.g. zq and 
z\, respectively. For that purpose we use the mapping, 

71(2) = iJ- — — 
V z - z 

In the next step we find a map that connects Z2 to the real 
axis by a semi circular arc. The inverse of this mapping, 
7^ 2 , brings Z2 to the real axis. 

where £ 2 = 7l (z 2 ) and a = |6l7^6 and b = |&|7»&- 
Iteratively, we apply this mapping to all the points 
23, . . . z n where in general 

£fc =7&_i °---°7i0fc)- 

In the third and last step we unfold the remaining part 
of the interior to the whole upper half-plane by the map 

with 

£n+l = 1U • • • 7i ( z o) 

The conformal map 'J' from the upper half-plane to 
the interior domain and from the lower half-plane to the 
exterior domain is then given by 

* = 7i _1 Ig 1 ' ' ' 7^ 7,7+1 ■ 

From the conformal map 'J we easily construct the 
map $ from the unit circle to the inclusion. In Fig. Owe 
illustrate how this is done. 
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FIG. 7: (7 yy evaluated along the positive imaginary axis (el- 
lipse) 



FIG. 6: Sketch of the construction of the conformal map from 
the exterior unit circle to the exterior of the inclusion. We 
choose a such that infinity is mapped to infinity. Similarly, 
for the interior map we choose a such that zero is mapped to 
zero. 



V. EXAMPLES 

In order to check the validity of our method we cal- 
culated the stress fields created by inhomogeneities with 
2 different geometries: an ellipse with semi axes 0.9 and 
vl + 0.9 2 (aspect ratio of about 1.5) and a smoothed 
triangular curve 1.752 + ^ (see fig. 2]). In the case of 
the elliptical inhomogeneity we compared our method to 
the known analytical solution which was first obtained by 
Hardiman [l[ in 1954. In the example below, the bound- 
ary conditions at infinity were set to er^ = -lAt. The 
shear moduli used were [i\ = 1-^ for the inhomogeneity 
and ii2 = for the matrix. The Poisson ratio was 

taken to be v = 1 for both inhomogeneity and matrix. 
In figs. [5] and [7] we can see the components of the stress 
field calculated outside the ellipse. The blue line is the 
stress calculated using Hardiman's solution and the red 
spots corresponds to the values obtained by our method. 
Similarly, we have calculated the stress field outside the 
triangular-like inhomogeneity (figs [8] and [9|) . 



VI. CONCLUDING REMARKS 

In comparing our approach to other available algo- 
rithms, for example finite elements approximations to 
the equations of linear elasticity, we should stress that 
our approach works equally well for compressible and in- 
compressible materials, There is no problem in taking 
the incompressible limit as the Poisson ratio approaches 
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FIG. 8: a xx evaluated along the positive real axis (triangle) 
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FIG. 9: dyy evaluated along the positive imaginary axis (tri- 
angle) 
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1/2. This is not the case for finite elements methods. 
While the examples shown above worked out very well, 
indicating that the proposed algorithm is both elegant 
and numerically feasible, unfortunately it deteriorates 
very quickly when the shape of the inhomogeneity de- 
viates strongly from circular symmetry. The difficulty in 
matching the two conformal maps is significant, as can 
be gleaned from from Figs. 2 and 3. One could think 
that the problem could be overcome in principle by in- 
creasing the numerical accuracy, but in practice, when 
the inhomogeneity has horns, spikes or deep fjords, the 
difficulties becomes insurmountable. Similar difficulties 
in another guise are however expected when any other 



analytic or semi-analytic method is used, leaving very 
contorted inhomogeneities as a remaining challenge for 
elasticity theory. 
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APPENDIX A: EXPLICIT CALCULATION 

Starting from ([36|) and (j37|) we write the series expansions in a compact form: 

OO DO a / m OO 

^(0 = ^(0 = =£**<*. (ai) 

fe=0 k=l k=-oo 

1 1 i / v OO 

= e = E ^- fc > == = E c ^ k . (A2) 

k— — oo k— — oo v ' — oo 

were we have eliminated the zero terms in ip 1 because of the gauge freedom. Substituting: 



E^c fc + (E ^oE^c^+E^r 



—k 

Vk 1 

k=0 n= — co k=0 k=l 







= J2 ^ k + ( E c «-") E + E ^ k + <&> + + ^7 ■ 

k— — oo n= — oo k— — oo fc= — oo ^ ' 

We can also expand: 

^ w + + — = 

CToo £H p ^oo (Too ^1 _ V f() n /AO) 

V*V n= — oo 

Substituting and changing the indices of summation: m = 1 — k + n which leads to n = m + k — 1 we get : 

oo oo oo — 1 

£^c+£ E kb m+k ^i\ m + e ^iu m = 

m— A;— 1 771= — oo m— — oo 



E E E kc m+k - 1(p fcj m + e ^- m ^ m + E /»^ m - ( A4 ) 

m— — oo k— — oo m=— oo m— m— — oo 
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In order to find the relations between the coefficients we need that both sides will be expressed in the same Fourier 
harmonics. Therefore, we need to expand: 

w (C) = x(A(0). 

In a Fourier series. The condition for this to be expanded in a Fourier series is that ui(Q is L 2 (i.e. \u(e ie )\ 2 d9 < oo 
on the segment [— If that is the case, we can expand: 

oo 

w(o - E a ^ k ■ 

k— — oo 

Actually, it is found to be more convenient to expand powers of ( in a Fourier series: 

oo 

^(cr = E a ^ n > ( A5 ) 

n— — oo 

were ( — e l6 and \cu(e l6 ) m \ 2 dO < oo on the segment [— 7r, tt]. 

oo oo oo —1 

X>»c n+ E E kb n+k w k *c+ E ^-nC = 

n—0 k—1 n— — oo n— — oo 



= E t E ^m a «>™ + E fc ( E C m+fc-l»n,m Vfc* + E V>-m a n,m + E ./™ a ™,™K" ■ (A6) 



n— — oo m— — oo k— — oo m— — oo m— 

Define: 



and 



Substituting: 



-Bn.fc — ^ ' Cro+fc— l a n,m (A7) 



m— — oo 



Cn — E fm a n,m (A8) 



m— — oo 



oo oo 



E^«c" + E E fcw-i^*c"+ E ^-nC" = 



n—0 fc— 1 n— — co 



oo —1 —1 oo 

- E t E ^ a «.™ + E kB ^i* + E + c nK n • (A9) 

n— — oo m= — oo k— — oo m— 

Using the linear independence of the £™ 's with respect to the Fourier integral, and 'cutting' the infinite series at some 
number N, we get a set of linear equations which is of the form: 

Miv = c (A10) 

Where v is the vector of coefficients (i/j's and ip's), M is a (47V + 2) x (2N + 1) matrix of constants and c is the C„'s. 
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Next, we substitute the expansions in the continuity equation for the displacement: 
Mi A'(C) M2 

Substituting: 

- oc oo oo oo 

± [Kl j2ric k -( E &nC)E^c 1 - fc -E^*c- fc ] = 



Ml 

/c— n— — oo fc— 1 fc— 1 



— K 2 2^ ^ W ~ 2^ ) fc ^ U - 1^ U + «Wl^ - ===<£i . (All) 

fc= — oo n= — oo fc= — oo fc= — oo > ' 



We can also expand: 



2 *M 2 , 



(Too $(w) CToo (Too -Fl „ ,, 10 , 

«2-Fi '■' - - F i = > QnW ■ (Al2) 



n— — oo 



Substituting and changing the indices of summation: m = 1 — k 4- n which leads to n = m + fc - 1: 

j OO CO oo oo 

-[«iX>£C fc -£ E fc& ro+ ,-i^C m -E^*C- fc ] = 



' fe=0 fe=l n=-oo fe=l 



-1 



= -[ K2 e E E fow*-^"™ - E ^-> m + E ^ m ]- (Ais) 

M2 — — — — ' 

m— — oo k— — oo 771= — co m— m=— oo 

Expanding cj(C) &s before and defining: 

oo 

-Dn = E 9ma n , m , (A14) 

m— — oo 

we get: 

- oo oo oo oo 

-kE^c"- E E fc6 »+*-^*£"-E^"] = 

n— n— — oo fe=l n— 1 

^ oc —1 —1 oo 

= E t K2 E rfn a n,m - E kB ^,W*k ~ E ^- m a «,m + Aj]C™ • (A15) 

^ 7i — — oo m=— oo fc— — oo m— 

When 'cutting' the infinite series in the same way, we get again a matrix equation (with the same dimensions) of the 
form: 

M 2 v = d (A16) 
Combining equations (|A10[) and (|A17|) we get a (AN + 2) x (47V + 1) matrix equation: 

Mv = e (A17) 

where 

M = Mi © M 2 (A18) 
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and 



cffid 
I 



(A19) 
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